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Simulations of the two-dimensional Ising and 3-state Potts models at their critical points are 
performed using the invaded cluster (IC) algorithm. It is argued that observables measured on a 
sub-lattice of size I should exhibit a crossover to Swendsen-Wang (SW) behavior for I sufficiently 
less than the lattice size L, and a scaling form is proposed to describe the crossover phenomenon. 
It is found that the energy autocorrelation time t e (1,L) for an / x I sub-lattice attains a maximum 
in the crossover region, and a dynamic exponent z ic for the IC algorithm is defined according to 



T E ,max ~ L z . Simulation results for the 3-state model yield z ic = .346 ± .002 which is smaller 
than values of the dynamic exponent found for the SW and Wolff algorithms and also less than the 
Li-Sokal bound. The results are less conclusive for the Ising model, but it appears that z IC < .21 
and possibly that r e . max ~ log L so that z IC = — similar to previous results for the SW and Wolff 
algorithms. 

I. INTRODUCTION 

Monte Carlo (MC) methods used to simulate classical spin systems, such as Potts models, fall primarily into 
two broad classes: local-update algorithms and cluster algorithms. Algorithms with local update rules, such as the 
Metropolis algorithm, provide an efficient means of simulating these spin systems in non-critical regions. Near a 
second-order phase transition, however, where long-range correlations are present, relaxation times increase rapidly 
with system size. This phenomenon, known as critical slowing down, may be characterized by a dynamic exponent 
z according to r ~ L z , where r is the autocorrelation time measured at criticality (roughly, the time necessary to 
generate a statistically independent configuration) and L is the system size. Local-update algorithms typically have 
values of z slightly greater than 2 and, therefore, are impractical for simulating large systems near a critical point. 

Cluster algorithms, on the other hand, such as the Swendsen-Wang (SW) algorithm, employ non-local update 
moves, flipping clusters of spins of linear extent comparable to the correlation length. This technique significantly 
reduces critical slowing down, and thus makes cluster algorithms preferable for simulating spin systems near a critical 
phase transition. Recent numerical estimates of the SW dynamic exponent for two-dimensional ferromagnetic q-state 
Potts models are z sw w .25 § for the Ising (q = 2) case and z sw « .52 j| for q = 3. 

Although there currently exists no theoretical means by which the dynamic exponent of a SW-type algorithm may 
be calculated, there is a rigorous lower bound. The Li-Sokal bound, as it has come to be known, states that 
z sw > a/u for q-state Potts models, where a and v are the usual static critical exponents for the specific heat and 
correlation length, respectively. We note that the numerical values given above are consistent with this bound, since, 
in two dimensions, ajv = (log) for q = 2, and a/v = 2/5 for q — 3. For the remainder of this paper we will continue 
to focus on the Ising and 3-state Potts models in two-dimensions, since these are the two most carefully studied cases. 

We now turn to the invaded cluster (IC) ]7|-|i"(| algorithm, a recent approach based on invasion percolation, for 
simulating equilibrium critical points. This algorithm has the unique property that it "self-organizes" to the critical 
point. Therefore, no a priori knowledge of the critical temperature is required; instead, T c is an output of the 
algorithm. In addition, due to an intrinsic negative-feedback mechanism, the IC algorithm equilibrates very quickly 
in the sense that thermodynamic quantities are measured to be near their equilibrium values within a few MC steps 
(after starting, say, from a completely ordered state). 

Initial studies [H] seemed to indicate that the IC algorithm suffers no critical slowing down for the Ising model. For 
both d = 2 and d — 3, the integrated autocorrelation time t £ was observed to decrease with L, while r m remained 
constant (within error bars), where e is the energy per spin, and m is the fraction of spins in the largest cluster. 
(We omit the usual "int" subscript on r, since we deal almost exclusively with integrated as opposed to exponential 
autocorrelation times.) The decrease of r e with L was also observed for 3- and 4-state Potts models in two dimensions, 
but in these cases critical slowing was evident in the behavior of r m . In Ref . [pj the dynamic exponents were estimated 
to be z m « .28 for q = 3 and z m w .63 for q = 4, each of which is less than the Li-Sokal bound on z sw for its respective 
value of q. 
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In this paper, we repeat these studies for q = 2 and q = 3 in two dimensions using larger lattice sizes and an improved 
method ||] of estimating r. We also investigate the L-dependence of the "specific- heat-like" quantity c(L) = L d var(e). 
(In the canonical ensemble c is the specific heat, but the same is not true for the IC ensemble.) In addition we measure 
r e (l,L), the integrated autocorrelation time for the energy per spin s(l,L) measured on an I x I sub-lattice of the 
whole L x L lattice, as well as c(l,L) = l d v&r(e(l, L)). While the motivations for these experiments are discussed in 
more detail in Sec. [jlj, the central idea is that we expect the negative feedback to diminish for length scales I < L, 
leading to sub-system behavior that differs from that of the whole system. We argue in Sec. [II that, for I sufficiently 
less than L, we should observe a crossover to SW behavior for observables measured on a sub-lattice of size I. Upon 
investigating the crossover region, we find that c(l,L) has a simple scaling form and give an estimate of the length 
scale at which c(l,L) crosses over to its SW analog, namely, the specific heat for an / x / system. 

The results are less conclusive in the case of the dynamic variable r e (l,L), but it appears that crossover to SW 
behavior does occur and that the crossover length for r e (l,L) differs from that for c(l,L) and is likely the same for 
q = 2 and q = 3. In addition, the crossover phenomenon leads to a maximum in r e (l, L) as I is varied for a given L. 
We argue in Sec. Ill that the dynamic exponent z IC of the IC algorithm is appropriately defined by r e , max ~ L z ' . 
For the 3-state Potts model we give a numerical estimate for z lc that is smaller than a recent estimate || of z sw 
and also less than the Li-Sokal bound on z sw . For the Ising model the results are less conclusive, but it appears that 
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< .21 and possibly that r e 



log L so that 



— similar to the state of affairs for the SW algorithm 



The remainder of this paper is organized as follows. In Sec. [n] we pr ovide some background on the invaded 
cluster algorithm and discuss results of previous IC simulations. In Sec. Ill we discuss the crossover phenomenon in 
greater detail and propose a scaling form to relate observables measured on sub-lattices using the IC algorithm to 
corresponding quantities for the SW algorithm. We describe our IC simulations of the Ising and 3-state Potts models 
in Sec. IV and discuss the results in Sec. M. Sec. VI contains our conclusions. 



II. INVADED CLUSTER ALGORITHM FOR CRITICAL POTTS MODELS 

In order to understand how the IC algorithm works, we first review the SW algorithm for Potts models. The 
ferromagnetic g-state Potts hamiltonian is 

(id) 

where <7j € {0, 1, . . . , q — 1} and the sum is over nearest-neighbor spin pairs. (Note that the Ising model is just the 
special case q = 2.) 

Given an initial spin configuration, the SW algorithm proceeds as follows: First, satisfied bonds are occupied with 
probability p = 1 — e~^, with (3 = 1/T, where a bond joining two spins i and j is defined to be satisfied if and only if 
Cj = aj. Unsatisfied bonds are never occupied. Next, clusters of spins connected by occupied bonds are identified, and 
each cluster is "flipped" , i.e., independently and uniformly assigned a new random spin value from {0, 1, . . . , q — 1}. 
(Note that a cluster can consist of a single spin.) Finally, after statistics have been collected, occupied bonds are 
erased and the whole process is repeated. It can easily be shown that the SW algorithm satisfies detailed balance for 
the canonical ensemble. 

The IC algorithm uses invasion percolation to generate the spin clusters to be flipped. Given an initial spin 
configuration, the first step is to assign a random order to the bonds of the lattice. The bonds are then examined, and 
satisfied bonds occupied one at a time in this order. If a bond joining two clusters is occupied, they are combined into 
one. Cluster growth continues until some stopping condition is fulfilled. In this paper, we consider the topological 
spanning condition, which dictates that growth be stopped as soon as some cluster winds around the system in one 
of the d directions. As soon as spanning is detected, clusters (including the spanning cluster) are flipped exactly as 
in the SW algorithm, statistics are collected, bonds erased, and the process repeated. 

To understand why the IC algorithm self-organizes to the critical point, we define / to be the ratio of the number of 
occupied bonds to the number of satisfied bonds when some cluster first spans the system. It has been argued @,||Jl0) 
that, as the system size L approaches infinity, the distribution of / approaches a delta function at p c = 1 — e _/3c , where 
[3 C is the inverse critical temperature. Though not a rigorous proof, the argument proceeds as follows. First, we note 
that p c is the threshold for percolation on the satisfied bonds of a critical spin configuration jTlfl . Thus, given a spin 
configuration that is typical of the critical point, the fraction / of satisfied bonds that must be occupied to achieve 
spanning is close to p c . Second, we observe that each iteration of the IC algorithm is identical to an iteration of the 
SW algorithm with p = f. Therefore, performing an iteration of the IC algorithm on a critical spin configuration is 
equivalent to performing an iteration of the SW algorithm with p p c , and thus the system will remain near the 
critical point. 
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If, instead, the system is started in the low-temperature phase, the number of satisfied bonds will be larger than 
is typical of T c . Therefore, a smaller fraction / will need to be occupied to achieve spanning. In this case, an IC 
iteration is equivalent to a SW iteration with p < p c , i.e., T > T c , and therefore the system is pushed toward T c from 
below. Similarly, if the system is started in the high-temperature phase, it is pushed toward T c from above. Thus, in 
summary, a system in a non-critical state is pushed toward criticality, while a system in a critical state remains near 
criticality with T fluctuating about T c . 

Because of this negative-feedback mechanism, the IC algorithm self-organizes to the critical point with no a priori 
knowledge of T c . Instead, T c is obtained as an output of the algorithm, via the relation (/} = 1 — e~ x l Tc . For example, 
results of IC simulations for the 2d Ising model yield an estimate of T c w 1.1355 when extrapolated to L = oo, as 
compared to the exact result for an infinite system, T c — 1.1346 . . .. 

Since every configuration generated by the IC algorithm must contain a spanning cluster, it is clear that the 
algorithm does not sample the canonical ensemble for a system of finite volume. We refer to the stationary distribution 
sampled by the IC algorithm as the IC ensemble. If we assume that as L — > oo the distribution of / approaches a delta 
function at p c and that the volume fraction of the spanning cluster goes to zero, local observables such as internal 
energy and magnetization will approach their infinite- volume critical values in this limit. Simulation results RPIJIol 
support this hypothesis. For example, results of IC simulations for the 2d Ising model yield an estimate for the energy 
per spin of e c ~ —1.706 JtJ when extrapolated to L = oo, as compared to the exact result for an infinite system, 
e c = -1.7071.... 

Some finite- volume fluctuations in the IC ensemble, however, are very different from those in the canonical ensemble. 
For example, in the canonical ensemble the quantity c(L) = L d var(e) is the specific heat which diverges as L a l v at 
the critical point — a logarithmic divergence for the Ising model in two dimensions. In the IC ensemble, however, 
c(L) is observed Q to increase roughly linearly with L for the 2d Ising model. These differences can be traced to 
fluctuations in the effective temperature (measured by /) in the equilibrium state. In the next section, we will examine 
more closely the roles played by temperature fluctuations and the negative-feedback mechanism in determining the 
properties of the IC algorithm. 

III. CROSSOVER TO SWENDSEN-WANG BEHAVIOR 

As described in the previous section, the negative-feedback mechanism, inherent in the IC algorithm, drives the 
system to criticality by effectively adjusting the temperature after each iteration. As previously mentioned, this 
mechanism leads to differences in the L-dependence of several dynamic and static quantities from that observed for 
the SW algorithm. Now, however, we consider an Z x Z sub-lattice within the L x L lattice. Since the energy of a 
sub-system of size Z <C L is weakly correlated with that of the whole system, the negative feedback mechanism is less 
effective for the sub-system. A "warm" (relative to T c ) sub-system in a "cool" system will be further warmed by the 
next IC iteration. As a result, the energy autocorrelation time for the sub-system may be longer than for the whole 
system. 

The observation that an iteration of the IC algorithm is equivalent to an iteration of the SW algorithm with p = f 
provides further insight into IC dynamics. In particular, if the distribution of / approaches a delta function as L — > oo, 
then any finite sub-system of an infinite system will behave exactly as it would under SW dynamics. In short, for I 
sufficiently smaller than L, the sub-system doesn't "know" it is being updated by the IC and not the SW algorithm. 
Thus we expect that, in the limit L — > oo, all static and dynamic quantities measured on a sub-system of finite size 
I will approach the values measured for the SW algorithm for a sub-system of the same size. It also follows that, 
for fixed L, there is a crossover from SW to IC behavior at intermediate values of I. For example, the integrated 

SW 

autocorrelation time t e (1,L) for the energy per spin in a sub-system of size Z should initially increase with I as l z 
for I <C L, reach a maximum, and then decrease as I is increased further into the range where the negative- feedback 
mechanism becomes significant. 

We can estimate the length scale at which the crossover occurs as follows. In the canonical ensemble, the temperature 
uncertainty of the critical region scales with system size L as ST ~ L~ x l v . In the IC ensemble, however, temperature 
fluctuations are governed by the negative-feedback mechanism as described above. In Ref. jl(J it was found that 
the standard deviation of / scales as <r(f) ~ L~ h with b rs .46(.30) for q — 2(3) as compared with \/v = 1(6/5) 
respectively. Since the SW algorithm samples from the canonical ensemble and since we expect SW behavior for sub- 
systems of size Z -c i, crossover between the IC and SW regimes should occur when the temperature uncertainties 
from the two sources are comparable. Thus, for a given L, we expect crossover at a sub-system size Z c given by 

Therefore, in light of the previous arguments, we hypothesize that the crossover from IC to SW behavior may be 
described by the scaling relationship 
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A IC (l,L) = A™(l)F A (l/Ly), (3.1) 

where A^(l) is any observable measured for the SW algorithm on an / x I sub-lattice immersed in an infinite system, 
A lc (l,L) is the same observable measured for the IC algorithm on an I x I sub-lattice of an L x L lattice, F A is a 
scaling function with the property that F A (x) — ► 1 as x — > 0, and y — bv. A word of caution concerning boundary 
conditions is in order here. We define Aq W (L) to be an observable measured for the SW algorithm on an L x L 
lattice with periodic boundary conditions. Although we expect A™ (I) / Aq W (I) to approach a constant for I — > oo, 
the constant will in general not be exactly 1 due to the different boundary conditions. 

In this paper, we also seek to define a meaningful dynamic exponent z lc for the IC algorithm that may be compared 
with exponents for other algorithms as well as with the Li-Sokal bound. It is not obvious how to do this since the 
energy autocorrelation time r e for the entire system was observed || to decrease with L. However, this is not the 
whole story since we have argued above that correlations between successive IC configurations should decay more 
slowly on length scales I < L. Since we would like z lc to describe the L-dependence of the slowest mode, we suggest 
that it is the maximum value of T e (l,L) attained for a given L that is relevant. Thus we define z lc according to 

IC 

T e,max ~ L Z . (3-2) 

In the next section, we describe simulations designed to measure z lc and to test the scaling hypothesis (Eq. ( |3.l| )) for 
the static variable c(l, L) as well as for the dynamic variable t s (1, L). 



IV. DESCRIPTION OF SIMULATIONS 



We used the invaded cluster algorithm with the topological spanning rule to simulate the two-dimensional Ising 
and 3-state Potts models at their critical points for systems ranging in size from L = 32 to L — 1024. Starting from 
a completely ordered state, we performed a number of relaxation steps to allow the system to reach equilibrium at 
T c and then collected data for four observables: the energy per spin e, the ratio / of occupied to satisfied bonds, 
the fraction m of spins in the largest cluster, and the susceptibility x, given by the sum of the squared cluster sizes 
divided by the total number of spins. In addition to measuring the mean value and variance for each observable, we 
also measured the autocorrelation function and used this to calculate integrated autocorrelation times. For a given 
observable A, the (normalized) autocorrelation function at a given time step t can be calculated from a sequence of 
n MC measurements {A(j),j = 1, . . . ,n} according to 

r (f) _ EU(A(j)-(A))(A(j+t)-(A)) 

TA{t) = YUM) -(A)? ' (41) 

where (A) is the mean value of A. 

The integrated autocorrelation time for the observable A is defined by 



TA -2 



1 OO 



Obviously, in practice, the sum must be truncated at some reasonable value of t. Following the recommendation of 
Ref. we define 

T A (t A )^- + J2r A (t) (4.3) 



4=1 



and choose the cutoff t* A to be the smallest integer such that t* A > KT A (t A ), where k is a constant whose value 
will be discussed shortly. If the autocorrelation function has the scaling form T A (t) = G(t/T exp ), where T oxp is the 
exponential autocorrelation time, then choosing the cutoff in this manner will insure that T A {t* A ) is proportional to 
t a . Thus estimates of z A will not be biased by truncating the sum at t = t* A , 

One also would like the values of T A {t* A ) to approximate t a as accurately and precisely as possible, and here there is 
a tradeoff between excluding noise and including as much of the signal as possible. In Ref. [|| it is shown that if T A (t) 
is roughly a single exponential, then choosing a value of k in the range 4-6 would achieve the optimal compromise for 
n/r in the range 10 4 -10 6 that we used in our simulations. However, although T A (t) is well approximated by a single 
exponential in the case of the SW algorithm, this is not true for the IC algorithm H. For this reason, and since we 
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are willing to accept slightly larger statistical uncertainties in order to reduce systematic errors, we used k = 10 in all 
our calculations. 



In light of the discussion in Sec. Ill, we also collected data for several sub-system sizes for each L. Here we 
concentrated on the energy per spin, measuring (e(l,L)), var(e(Z,L)), and r e (l,L) for sub-system sizes ranging from 
/ = 1 to I = L/2 for each L. The sub-systems are squares all sharing a single corner of the lattice. 

Error bars on all quantities were calculated using the blocking method. Each run was partitioned into k contiguous 
blocks of n MC steps each, and the individual blocks treated as independent runs. Although this is an approximation, 
it will be a good one provided that n is large compared to the system's longest relaxation time. As an example of 

(i) 

the blocking method, we obtain the value for, say, r m by first calculating r„ for each block i. We then calculate the 
mean 



fc 

h ^ 



^ (4.4) 

K, 

i=l 

and its standard error 



a(Tm) = V k(t—T) (45) 

and report the result T m — Tm ± &(jm)- 

For both q = 2 and q = 3, one long run was initially performed for each lattice size L, and the blocking method 
implemented as just described. (Three independent runs were performed for the case q = 2, L = 1024.) The number 
of blocks used ranged from 100 for L — 32 down to 10 for L = 1024, and the number n of MC steps per block ranged 
from 5 x 10 3 to 1 x 10 5 . In each case, n was greater than the longest observed autocorrelation time for the given system 
by at least a factor of 10 3 (10 4 for the smaller lattices), thereby making the assumption of independent blocks, used 
in calculating the error bars, a reasonably good one. The number of equilibration steps performed at the beginning 
of each run also exceeded the longest observed r by a factor of 10 3 in all cases. 

In these initial runs, we collected data for sub-system sizes I € {1,2,4,..., L/2} as well as for the whole system 



(I = L). In order to estimate z ic as defined in Eq. (3.2), we sought to obtain an accurate value for r £imax (L) for 



each L. Therefore, once we had learned, from the initial runs, the approximate sub-system size Z ma x a t which r e (l,L) 
attains a maximum, we then performed between one and three additional independent runs for each system and 
collected data for evenly spaced values of I near our rough estimate of Z max . The entire experiment required about 5 
months of CPU time on a single processor of a dual-processor 266MHz Pentium II Linux workstation. We used the 
machine-supplied random number generator coupled with a shuffling procedure as described in Ref. [L0| . 



V. DISCUSSION OF RESULTS 

The first quantity we consider is the static variable c(L) = L d va,r(e) (see Table Q). In the canonical ensemble, c(L) 
is the specific heat which diverges, at the critical point, as logL for q = 2 and as L a l v with ot/v = 2/5 for q — 3. We 
see from Fig. ffl, however, that the situation is quite different for the IC ensemble, as first observed in Ref. 0. We 
assume that the asymptotic behavior is given by a power law c(L) ~ L w and fit a line to a plot of log 10 c{L) vs log 10 L 
as shown in Fig. 0. Note that where error bars are not visible they are smaller than the symbol height. As is always 
the case when trying to ascertain asymptotic behavior from simulations at finite L, there can be some debate as to 
which, if any, data points should be omitted from the fit because of corrections to scaling. Here and elsewhere in our 
analysis, we proceed by dropping points one at a time in order of increasing L until either 1) a reasonably good fit is 
obtained, 2) the fit ceases to improve significantly with further cuts, or 3) we are left with only three data points. We 
employ standard, weighted x 2 fitting, using the confidence level [] (CL) as our goodness-of-fit measure, and consider 
a fit to be "reasonably good" if CL > 10%. 

For the Ising model, a fit to the last four data points (L > 128) yields w — 1.020 ± .003 (CL = 16%) in agreement 
with the observation w rj 1 reported in Ref. M . In the case of the 3-state Potts model, a fit to the last three points 



x The confidence level is the probability that a \ 2 as poor as the measured value would occur, assuming that the underlying 
model is correct and that the measurement errors are normally distributed ]l| . (Note that confidence level is denoted by the 
symbol Q in Ref. @.) 
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gives w = 1.313 ± .008, but, because of the poor confidence level (2%) and the upward curvature visible in the data, 
this value should probably just be regarded as a lower bound on w for q = 3. We emphasize that the error bars 
on these and subsequent exponent estimates are purely statistical in nature and do not reflect the uncertainty of 
extrapolating to infinite system size. 

Turning now to the dynamic variables (see Table |nj) , we plot the logarithms of the autocorrelation times r E , ry, 
T m , and t x vs. log 10 L for q = 2 in Fig. || and for q = 3 in Fig. ||. For q = 2 we find that t e and tj decrease with 
L (perhaps in a rather complicated fashion) and r m and r x remain approximately constant for the range of L values 
used in our simulations. These results are in agreement with initial observations Q that led to speculation of no 
critical slowing. For q — 3, however, we observe critical slowing in the behavior of r m and r x . Fits to the data for 
L > 64 yield z m = .191 ± .004 (CL = 25%) and z x = .206 ± .005 (CL = 20%), but it seems likely that r m - t x in the 
asymptotic limit. We note that the value of z m is somewhat smaller than the previous estimate z m re .28 

Next we consider T e (l,L), the integrated autocorrelation time for the energy per spin e(l,L) measured on an I x I 
sub-lattice of the whole L x L lattice. As expected from the discussion in Sec. II], we see that as I is increased for 
a given L, T e (l, L) increases, reaches a maximum, and then decreases, as shown in Fig. ^ for q — 2 and in Fig. [| for 
</ ■>■ 

To find r e max and its location / max for a given L (see Table HI), we fit a parabola to the region of the curve near 
the maximum. In order to do this objectively, we began, for each L, by omitting the data point with the smallest 
value of r e and performing the fit. We then dropped the point with the next smallest r E , re- fit, and continued in 
this fashion until 1) a CL of 50% or greater was obtained and 2) the values of T e max and l max that were obtained by 
dropping an additional point remained within error bars of the current best-fit va lues . 

We then attempted to determine the dynamic exponent z lc , defined in Eq. (3.2), by fitting a line to a plot of 
log 10 T E)max (L) vs. log 10 L for q = 2 and q = 3 . The results are shown in Fig. [6] along with results for rf w taken from 
Baillie and Coddington S for q = 2 and from Salas and Sokal j| for q = 3. We note that the increase of T S;laax (L) 
with L for q — 2 is the first observation of critical slowing for the IC algorithm in the case of the Ising model. 

For q = 2 the autocorrelation times for the IC algorithm are nearly the same as for the Wolff algorithm and smaller 
than those for the SW algorithm by a factor of about 2 (see Table |lj]), but the L-dependence is similarly obscure. A 
good fit to a power law could not be obtained for the IC data for q = 2. The line shown in the figure, having slope 
re .21 is the best fit for L > 64, but it clearly does not describe the data very well. The best power-law fit to the 
SW data for 64 < L < 512 (also shown in the figure) yields z IC re .25 as reported in Ref. Although the fit is 
considerably better than that for the IC data, the CL is still poor (< 0.1%), and the better fit might be primarily 
due to the absence of data for L = 1024 in the SW case. 

Since it has been suggested |l3|] that r^ w increases logarithmically with L rather than as a power of L, Baillie and 
Coddington also fit their data to a logarithm with somewhat better results (CL = 13%). For the IC algorithm, a 
logarithmic fit is still atrocious, albeit somewhat better than the power law. Later in this section we consider the 
possibility that we have underestimated the error bars on r 6imax , which, of course, could result in a poor fit even if 
the underlying model had been correctly identified. Still, even the general trend in the data is difficult to discern, 
indicating that corrections to scaling are probably significant for the system sizes studied here. Therefore, we conclude 
that high precision data for larger lattices are needed before a more definitive statement can be made concerning the 
asymptotic behavior of T £iinax . 

For the 3-state model, however, the picture appears to be somewhat clearer. Although a slight downward curvature 
in the data is visible in Fig. ^, a good fit (CL = 69%) to a power law is obtained for 256 < L < 1024, yielding 
z lc = .346 ± .002. This is to be compared with the value of zf w = .515 ± .006 obtained by Salas and Sokal for the 
SW algorithm with 128 < L < 1024 (CL = 80%). For the Wolff algorithm zf olft = .57 ± .01 was reported in Ref. §. 
While there is no guarantee that t E:Wl3X is the system's longest relaxation time, it is interesting that z lc is significantly 
smaller than zf w and also less than the Li-Sokal bound (z sw > a/v = 2/5). 



Now we proceed to test the scaling hypothesis (Eq. fl3.1| )) presented in Sec. [II. In order to do this, we first need 
to determine the exponent b defined by a(f) ~ L~ b . We plot log 10 a(f) vs. log 10 L in Fig. ^ (the data are listed in 
Table |) for q = 2 and q = 3 along with the best-fit lines. A fit to all six q = 2 points yields b = .4781 ± .0006 
(CL = 54%), while, for q = 3 a good fit (CL = 68%) is obtiained for the last four points [L > 128), resulting in 
b = .3252 ± .0009. Thes e results are consistent with previous estimates M. 

We now test Eq. Q for the variable c IC (Z, L) = 2 d var(e(Z, L)). Since cf w (Z) - log / and we expect c^(/) ~ cg w (0 
(recall Cq W (Z) is the specific heat for an I x I lattice with periodic boundary conditions, and (I) is the specific heat 
for an I x I sub-lattice immersed in an infinite system) , we plot c IC {l,L)/ log 10 / vs. I / L y in Fig . ||, where y = b = .4781 
for the Ising model (v = 1). The observed data collapse provides strong support for Eq. ( |3.l| ). For q = 3, however, it 
was found in Ref. Q that the asymptotic form Cq W (1) ~ \ a l v doe s no t describe the SW data very well for the range 
of lattice sizes considered here. Therefore, we cannot expect Eq. (p3j) to provide a good description of the IC data if 
the asymptotic form is used for c™(l). Nevertheless, if we plot c v (l, L)/cq W (1) vs. l/L y , using the measured values 
of Cq W (Z) from Ref. (D and y = bv = .2711 for q = 3, data collapse is apparent in Fig. ^|, although perhaps a bit less 
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convincing than for the Ising case. 

We note that the curve in Fig. || extrapolates to about .85 on the vertical axis for l/L v = 0. As previously mentioned, 
we would expect this value to be 1 if SW data were collected for sub-systems immersed in larger systems so as to 
repro duce the boundary conditions applied to the IC sub-systems. Thus we conclude that our scaling hypothesis 
(Eq. (3.1)) does appear to be valid for the variable c (l,L) for q = 2 and q = 3. 

Although the static quantity c (I, L) seems to be well des cribe d by Eq. (3.1) with y = bv, the same is not true for 
the dynamic quantity r e (l,L). This is easy to see, since Eq. (3.1) predicts that the location Z max of the maximum in 
Figs. Q and || should scale as L v ; however, the plots of log 10 l max vs. log 10 L shown in Fig. [h] reveal that this is not 
the case — at least not if y = bv is required. For both q = 2 and q = 3 the slope of the best-fit line for 64 < L < 1024, 
shown in Fig. [l^, is approximately .62, although the confidence levels are poor. Unlike the situation encountered 
earlier in this section, when fitting r 6 lnax to a power law in L, the points seem to be scattered randomly about the 
best-fit line. Therefore, we suspect that / max does scale as a power of L, but that our error bars on Z max are somewhat 
underestimated. 

There are three aspects of our analysis that could lead to underestimates in the error bars on r £jmax and i max . 
First, there always exists the possibility that the assumption of normally distributed measurement errors is not valid. 
Second, the blocking method, used to calculate error bars on values of t e {1, L), treats successive blocks as if they were 
independent runs, an approximation that may not be entirely justified even though the block length was greater than 
10 3 r ei m ax in all cases. Finally, and probably most importantly, the error bars on r £ , max and i max resulting from the 
weighted \ 2 fit to a parabola are calculated by assuming that the measurements of r e (l,L) at different values of I for 
a given L are independent. This is clearly not a good approximation, since all of the sub-lattices extend outward from 
the same corner of the L x L lattice. Therefore, all the spins in a given I x I sub-system are also contained in every 
larger sub-system, and thus sub-systems for comparable values of I are highly correlated. 

In any case, it is still clear that Eq. ( |3.l| ) with y = bv cannot explain the data for t £ . Since the logic leading to 
the scaling form seems sound, we hypothesize that Eq. ( |3.1| ) does hold for r e but that the crossover length for r e is 
different from that for c so that y ^ bv in the case of r e . This seems plausible, since there is no a priori reason why 
the t herm odynamic argument by which we arrived at y = bv must apply to the dynamic quantity r e . Nevertheless, if 
Eq. (3.1) still holds for t £ , we can obtain the crossover exponent y from Fig. [l(] as described above. 

To test our scaling hypothesis for r s , we plot r e (l , L) / 't°q '{I) vs. l/L v with y = .6176 (y = .626) for q = 2 (q = 3) 
in Fig. |ll] (Fig. |l|). The values of if J (I) are taken from' Ref. J| for q = 2 and from Ref. S for q = 3. The data 
collapse is not terribly convincing in either case, but seems too good to completely rule out Eq. (3.1) as the correct 
asymptotic form. The fact that both curves extrapolate to about 1 for l/L y = provides further support for the 
scaling hypothesis. Still, it appears that additional tests are needed to confirm or disprove Eq. (3.1) for t £ . 



VI. CONCLUSION 



Using the invaded cluster (IC) algorithm with the topological spanning rule, we simulated the critical Ising and 
3-state Potts models in two dimensions for systems ranging in size from L = 32 to L — 1024. In accord with 
previous results we find that the ^-dependence of several static and dynamic quantities is very different from 

that observed for the Swendsen-Wang (SW) algorithm which samples from the canonical ensemble. In particular, the 
quantity c(L) = L d var(e) is not proportional to the specific heat and the integrated autocorrelation time r e for the 
energy per spin decreases with L. However, we find that the corresponding quantities t s (1,L) and c(l,L), measured 
for a sub-system of size I, exhibit a crossover to SW behavior for I sufficiently less than L. 

To describe the crossover phenomenon, we propose the scaling form A lc (l,L) — A™ (1)Fa{1/ L y ), where A™ (I) is 
any observable measured for the SW algorithm on an / x I sub- lattice immersed in an infinite system, A lc (l, L) is the 
same observable measured for the IC algorithm on an I x I sub-lattice of an L x L lattice, and Fa is a scaling function 
with the property that Fa(x) — > 1 as x — > 0. We have argued that the crossover exponent y should equal bv, where 
v is the usual correlation-length exponent and b is defined by u{f) ~ L~ b with / the ratio of occupied to satisfied 
bonds. We find that the proposed scaling form with y = bv provides a good description of our data for the static 
variable c(l, L), but is less successful for the dynamic variable T e (l, L), even if the possibility y ^ bv is admitted. 

In addition, we define the dynamic exponent z IC for the invaded cluster algorithm in terms of the maximum value 
T £ ,max attained for a given L according to r £ymax ~ L z . For q = 3 we find that z lc = .346 ± .002, which is smaller 
than recent numerical estimates of z sw and z Wolff and also less than the Li-Sokal bound on z sw . For q = 2 we also 
observe critical slowing, but the L-dependence of r £!max is less clear. It appears from our simulations that z lc < .21 
and possibly that r £jmax ~ logi (z IC = 0), but high precision data for larger lattices are needed before a more 
definitive statement can be made concerning the asymptotic behavior of T £ rnax for the Ising model. 
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TABLE I. c = L d var(e) and er(/) for the IC algorithm for the 2d Ising and 3-state Potts models, where e is the energy per 
spin, / is the ratio of occupied to satisfied bonds, and L is the lattice size. 



L c (q = 2) c (q = 3) a(f) (q = 2) a(f) (q = 3) 

32 3.288(3) 5.067(4) .05100(4) .06991(4) 

64 6.038(6) 11.239(7) .03657(3) .05508(2) 

128 11.87(2) 26.24(2) .02626(3) .04367(3) 

256 24.03(8) 63.0(1) .01886(4) .03482(4) 

512 48.4(3) 155.0(6) .01352(7) .02781(9) 

1024 99.7(6) 390(2) .00977(4) .0223(1) 



TABLE II. Integrated autocorrelation times for the IC algorithm for the 2d Ising and 3-state Potts models, e is the energy 
per spin, / is the ratio of occupied to satisfied bonds, m is the fraction of spins in the largest cluster, \ is the susceptibility, 
and L is the lattice size. 



q 


L 




Tf 




T x 


2 


32 


.546(1) 


.1828(7) 


.857(3) 


.798(2) 


2 


64 


.499(2) 


.1275(7) 


.852(3) 


.795(3) 


2 


128 


.443(2) 


.0845(7) 


.859(3) 


.802(3) 


2 


256 


.384(3) 


.070(3) 


.869(8) 


.807(7) 


2 


512 


.305(3) 


.065(3) 


.88(1) 


.82(1) 


2 


1024 


.260(3) 


.033(2) 


.90(1) 


.83(1) 


3 


32 


.832(2) 


.1835(5) 


1.303(4) 


1.208(4) 


3 


64 


.821(2) 


.1742(4) 


1.435(4) 


1.351(4) 


3 


128 


.753(2) 


.1369(5) 


1.629(7) 


1.552(7) 


3 


256 


.635(3) 


.1071(8) 


1.89(1) 


1.81(1) 


3 


512 


.516(5) 


.079(1) 


2.13(3) 


2.08(3) 


3 


1024 


.407(4) 


.053(1) 


2.38(5) 


2.32(5) 



TABLE III. Zmax and autocorrelation times for the 2d Ising and 3-state Potts models. L is the lattice size, Z max and tJ^^ are 
the location and height, respectively, of the maximum in Figs. | and |. r £ sw and r e Wolff are the integrated energy autocorrelation 
times for the SW and Wolff algorithms. 



q 


L 


/max 


' £ ,max 


SWa 


^Wolff b 


2 


32 


8.36(5) 


1.962(3) 


4.016(5) 


1.815(3) 


2 


64 


13.89(8) 


2.319(2) 


4.90(1) 


2.225(6) 


2 


128 


21.1(4) 


2.694(7) 


5.87(2) 


2.654(12) 


2 


256 


32.2(3) 


3.133(6) 


6.87(3) 


3.076(24) 


2 


512 


51.2(6) 


3.60(1) 


8.0(1) 




2 


1024 


75.9(9) 


3.90(2) 






3 


32 


9.53(7) 


4.206(8) 


13.28(6) 


8.76(4) 


3 


64 


15.9(2) 


5.55(1) 


19.5(1) 


13.08(16) 


3 


128 


25.1(2) 


7.17(1) 


28.5(1) 


19.5(3) 


3 


256 


37.0(3) 


9.003(7) 


40.8(2) 


27.7(8) 


3 


512 


63.2(9) 


11.46(5) 


58.5(6) 




3 


1024 


90(1) 


14.5(1) 


82.2(2) 




a From Ref. | 
b From Ref. | 


for q = 2 and Ref. | 


| for q = 3. 
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FIG. 1. log 10 c(L) vs. log 10 L for the IC algorithm, plotted for the 2d Ising and 3-state Potts models. Here, c(L) = i d var(e), 
where e is the energy per spin and L is the lattice size. The solid (dashed) line is a linear fit to the q = 2 (q = 3) data for 
128 < L < 1024 (256 < L < 1024) and has slope 1.020 (1.313). 
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FIG. 2. log 10 r vs. log 10 L for the IC algorithm, where r is the integrated autocorrelation time and 32 < L < 1024 is the 
lattice size, plotted for the 2d Ising model for the energy per spin e, the ratio / of occupied to satisfied bonds, the fraction m 
of spins in the largest cluster, and the susceptibility x- 
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FIG. 3. log 10 r vs. log 10 L for the IC algorithm, where r is the integrated autocorrelation time and 32 < L < 1024 is the 
lattice size, plotted for the 2d 3-state Potts model for the energy per spin e, the ratio / of occupied to satisfied bonds, the 
fraction m of spins in the largest cluster, and the susceptibility x- 
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FIG. 4. The integrated autocorrelation time t s (1,L) for the energy per spin e(l,L) measured on an i x ! sub-lattice of an 
L x L lattice, plotted vs. log 10 I for L = 128, 256, 1024 for the IC algorithm in the case of the 2d Ising model. 
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FIG. 5. The integrated autocorrelation time t s (1,L) for the energy per spin e(l,L) measured on an I x I sub-lattice of an 
L x L lattice, plotted vs. log 10 I for L = 128, 256, 1024 for the IC algorithm in the case of the 2d 3-state Potts model. 
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FIG. 6. log! 



, r e sw and log 10 - 



for the 2d Ising and 3-state Potts models. r e is the integrated energy autocorrelation 



time for the Swendsen-Wang algorithm on an L x L lattice and r, 
lines are linear fits to the data (see text for further details). 
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FIG. 7. log 10 cr(/) vs. log 10 L for the IC algorithm, plotted for the 2d Ising and 3-state Potts models. a(f) is the standard 
deviation in the ratio / of occupied to satisfied bonds and L is the lattice size. The solid (dashed) line is a linear fit to the 
q = 2 (q = 3) data for 32 < L < 1024 (128 < L < 1024) and has slope -.4781 (-.3252). 
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FIG. 8. c IC (7,£)/log 10 Z vs. l/L y for the IC algorithm, plotted for the 2d Ising model for 32 < L < 1024 and 8 < I < L/4. 
Here, c lc (l,L) = i d var(e(Z, L)), where e{l,L) is the energy per spin measured on an I x I sub-lattice of an L x L lattice and 
y = bv, where b = .4781 is minus the slope of the solid line in Fig. m and v = 1 is the correlation-length exponent. 



17 



3.0 



"i r 



"i r 



"i r 



2.5 - 



2.0 - 



o 



y 1.5 



o L=32 

• L=64 

v L=128 

t L=256 

□ L=512 

- L=1024 



1.0 - 



0.5 I I I I I I I I I 

5 10 15 20 25 30 35 40 

III] 



FIG. 9. c IC (7,L)/cg w (7) vs. l/L v for the IC algorithm, plotted for the 2d 3-state Potts model for 32 < L < 1024 and 
8 < I < L/4. Here, c lc {l,L) = Z d var(e(/, L)), where s(l,L) is the energy per spin measured on an I x 1 sub-lattice of an L x L 
lattice, Cq(1) is the specific heat for an I x I lattice, y = bu is the crossover exponent, b = .3252 is minus the slope of the 
dashed line in Fig. (?] and v = 5/6 is the correlation-length exponent. 
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FIG. 10. log 10 Zmax vs. log 10 L for the IC algorithm, plotted for the 2d Ising and 3-state Potts models. i max is the location 
of the maximum in Figs. ^| and || and L is the lattice size. The solid (dashed) line is a linear fit to the q = 2 (q = 3) data for 
64 < L < 1024 and has slope .6176 (.626). 
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FIG. 11. ri G a,L)/r e S o (0 vs. l/L v for the 2d Ising model for 32 < L < 1024 and 8 < I < L/4. r £ IC (/,L) is the quantity 
plotted in Fig. W, T^[l) is the integrated energy autocorrelation time for the Swendsen-Wang algorithm on an I x I lattice and 
y — .6176 is the slope of the solid line in Fig. fLOl 
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FIG. 12. t £ ic (/,L)/t £ s T(Z) vs. for the 2d 3-state Potts model for 32 < L < 1024 and 8 < I < L/4. t £ ic (Z,L) is the 



quantity plotted in Fig. ^, r^(l) is the integrated energy autocorrelation time for the Swendsen-Wang algorithm on an I x I 
lattice and y = .626 is the slope of the dashed line in Fig. [l^. 
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